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Abstract 



A fully massive Monte Carlo program to compute all four-fermion processes 
in e~^e~ collisions, including Higgs boson production, is presented. Leading 
higher order effects are discussed and included. 
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PROGRAM SUMMARY 



Title of program: NEXTCALIBUR. 

Program obtainable from: R. Pittau, Dipartimento di Fisica Teorica, Univer- 
sita di Torino, Via Giuria 1, 1-10125 Torino-Italy. 

Licensing provisions: none. 

Computer for which the program is designed and others on which it has been 
tested: DIGITAL-ALPHA (double and quadruple precision) and HP (double 
precision only). 

Operating system under which the program has been tested: UNIX. 
Programming language used: FORTRAN 90. 
Memory required: about SOOkbytes. 
number of bits per word: 32. 
Subprograms used: none. 

Number of lines in distributed program: 10274. 

Keywords: decaying vector-boson production, all massive four-fermion pro- 
cesses, electroweak background, initial state pt dependent QED radiation, 
recursive Dyson-Schwinger equations, multi-channel Monte Carlo approach. 

Nature of physical problem 

Four-fermion production is and will be investigated at e"'"e~ colliders in a 
wide range of energies. Fermion masses have to be included when studying 
Higgs boson production and processes with forward scattered electrons. A 
Monte Carlo program being able to take into account both signal and back- 
ground electroweak diagrams for four-fermion processes, without neglecting 
fermion masses, is therefore required. Higher order effects, such as the QED 
radiation evaluated at the proper scale, have also to be included, together 
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with the possibihty of generating photons with a non vanishing pt spectrum. 
For t-channel dominated processes, the correct running of uem must be im- 
plemented. 

Method of solution 

An event generator is the most suitable choice for a program to be able to 
deal with the above physical problems. For the matrix element evaluation, 
recursive Dyson-Schwinger equations that express the n-point Green's func- 
tions in terms of the l—,2—,...,{n— 1)— point functions are used. The 
Monte Carlo integration is performed by using an efficient self-optimizing 
multi-channel approach. 

Typical running time: about 250 weighted events per second, in double preci- 
sion on DIGITAL-ALPHA/21164 DS20 (500MHz), depending on the chosen 
physical process. 

Unusual features of the program: by changing compilation procedure, the 
program can run both in double and quadruple precision. 

LONG WRITE-UP 
1 Introduction 

Electron-positron collisions have provided an ideal testing ground for the 
Standard Model. In recent years LEP2 has been instrumental in this. Future 
colliders may take over this role. At these high energies, final states will be 
produced containing unstable particles, like a single VT-boson, a Higgs boson 
or pairs of vector bosons. The net result of this will be a four-fermion final 
state. It is then not surprising that event generators producing four fermions 
are in demand. The situation at the start of LEP2 has been described, 
amongst others, in the 1996 Yellow Report and the present status in the 
2000 LEP2 Yellow Report [0. 

It has been our aim to provide an event generator, which can produce all 
possible electroweak four-fermion final states in all possible configurations. 
Moreover, its accuracy should match the experimental requirements as much 
as is feasible. 
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In the early LEP2 days an event generator like EXCALIBUR fulfilled 
this purpose to a certain extent by producing efficiently all possible massless 
four-fermion final states. At present, the above aim can only be strived for 
by 

• including fermion masses, 

• improving the treatment of QED radiation, 

• taking into account scale-dependent corrections and higher order con- 
tributions related to unstable particles. 

In the first place, fermion masses are needed both for a description of all 
kinematical regions (e.g. forward scattering) and for the inclusion of Higgs 
particles. 

In the second place, the accuracy should be improved by extending the 
treatment of QED initial state radiation in two ways. A pt effect of emitted 
photons should be generated and a realistic scale for t-channel dominated 
processes should be built in. 

Thirdly, the accuracy of the generator can also be increased when large 
fermion- loop effects on vector boson propagators are taken into account. 

The incorporation of these three effects into one event generator led to 
a completely new program, NEXTCALIBUR. It is the purpose of this paper to 
describe the program and the built-in physics treatments. Section |] refers 
to the new recursive matrix element calculation, section ^ summarizes the 
phase space generation, whereas the following two sections describe the above 
mentioned accuracy improving physics effects. Then, sections are devoted 
to the actual program, i.e. its structure, compiling instructions, input and 
test run output. Finally, two appendices clarify the multiperipheral phase 
space generation and the adopted fermion-loop scheme. 

For those readers, who are interested in the specific physics results of 
NEXTCALIBUR and some comparisons to other evaluations we refer to ref. 0]. 

2 The Matrix Element evaluation 

The algorithm for the computation of the matrix element is based on the 
Dyson-Schwinger equations, a set of recursive equations that express the n- 
point Green's functions in terms of the l—,2—,...,{n — 1)— point functions. 
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For a detailed presentation of the algorithm as well as its FORTRAN imple- 
mentation, HELAC, we refer to ref. ||^. Here we would like to summarize the 
main advantages of it: 

1. Scattering amplitudes involving any Standard Model interaction can 
be treated, notably the full Higgs sector. 

2. Masses of particles are fully considered. 

3. Both unitary and Feynman gauges are fully implemented. 

4. Quadruple as well as arbitrary precision arithmetic is available, in order 
to study special kinematical configurations. 

One important point we would like to underline here is that the algorithm 
makes use of the chiral representation for fermions and therefore it is more or 
less equivalent in efficiency to the method used by EXCALIBUR for massless 
fermions. 

Because of its modular structure, the algorithm can easily incorporate 
higher order contributions to the scattering amplitude, and in fact several 
such contributions have been included. The user has to specify the corre- 
sponding flags and input parameters, as described is section ^ in order to 
study these higher order effects. 

3 The integration strategy 

The integration strategy in NEXTCALIBUR is the same adopted in ref. [Q], 
namely a multi-channel self-optimizing |Q approach. 

All possible peaking structures of the integrand are taken into account by 
the 19 different kinematical channels (plus momenta permutations) in fig. 1. 
The conventions to read off the peaking structure follows. Fermionic 

lines have an arrow, wavy lines represent photons while dashed lines are 
massive gauge bosons. Sohd lines connect topological equivalent points: a t- 
channel solid line means isotropic angular distribution between the connected 
fermions, while an s-channel line only gives rise to an s dependent behaviour, 
without affecting the peaking structure. 

The difference with the channels of ref. P] is that now all fermions are 
taken to be massive. 
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Figure 1: kinematical diagrams in NEXTCALIBUR. 
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Furthermore, the leading kinematical structures for Higgs boson produc- 
tion have been added. They are represented by the channels NONABl (Ml ,M2) 
and N0NAB4(M), where the masses can now also be the Higgs mass. 

Generally speaking, the modification from the old massless to the new 
massive channels required only trivial changes. An exception is represented 
by the channel MULTIl, for which we list, in appendix 0, the complete gen- 
eration algorithm. 

A last remark is in order. When dealing with final state electrons, strong 
numerical cancellations occur, that may degrade the accuracy of the phase 
space generation. In order to overcome this problem a number of tricks are 
applied all over the places in NEXTCALIBUR, all derived from the basic identity 



TTL^ 

qo-\q\ = 7 fr^ , where = ml, q= (go, q) ■ (1) 

[qo + W 

4 Treatment of the QED radiation 

In NEXTCALIBUR QED radiation is implemented via the Structure Function 
formalism, namely by convoluting the Born cross section together with QED 
Initial State Radiators 0. Two ISR photons are explicitly generated, by 
using pt dependent Structure Functions derived, at the first leading loga- 
rithmic order, for small values oipt @] • Our starting point is the convolution 

/'27r ]^ /'27r 

dxi dx2 dci dc2 — rf^i — / d(f)2^{xi, Ci)^{x2, C2) (To{s) , (2) 
/vr JO 2tt Jo 

where cxo is the Born four-fermion cross section, s the reduced center of 
mass energy of the event, after photon emission, ci^2 and 0i,2 the cosines 
of the directions (in the laboratory frame) and the azimuthal angles of the 
two emitted photons with respect to the incoming particles. Finally, the 
Structure Functions are given by 

$(x,c,g2) = 

exp {1/^(1-,.)}^ _ I 1 _ 



2^ ^ \l-c + 2^ ^ 1 + ' 2 ~ (1 + x)(l - c + 2^ 
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TT 
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The above formula can be obtained by unfolding the strictly coUinear Struc- 
ture Functions, with the replacement 



mi 

for the soft part, and 



1 



c + 2^ 



- 2- 



:i-c+2f )2 
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mi 



l+x 



1 — s 1 

4_£. 



:i + x)(l-c + 2^^2 



(5) 



(6) 



for the hard contributions Q With this choice, the whole expression gets 
proportional to — 1, after integrating over c, as it should be. The 

inclusive QED result is therefore recovered and, at the same time, the pattern 
of the photon radiation is exact for small values of pt- 

The generation algorithm is as follows. The initial state electrons, before 
radiating, have momenta given by 



p, = {E, PE, 0, 0),P2 = {E, -(3E, 0, 0) , = ^ , /3 = y 1 - . (7) 

Once Ci^2 are generated, together with the energy fractions Xi^2 and the az- 
imuthal angles 0i,2, the momenta of the two ISR photons are completely 
determined 



Qi = E(l -xi)(l,ci,sicos0i,sisin0i) , 

q2 = E{1 - X2)(l, -C2, S2COS02, S2sin</)2) 



l-c? 



(8) 



^Note that the replacement of eq. (g) is also performed in the higher order terms of 



eq. 
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and the invariant mass, after ISR, is given by 



(1 - Xi)(l - X2) 




[ciC2 - 1 - S1S2COS(01 - 02)] 



2 

P2-q2- (9) 

We then generate four-fermion events in the center-of-mass frame of i, where 
the new x direction after ISR is determined by the spatial part of the four- 
vector obtained by boosting ii in the center-of-mass frame. All generated 
four-fermion momenta are subsequently boosted back to the Laboratory 
frame. 

In NEXTCALIBUR there is also the possibility of generating strictly collinear 
radiation. In this case the relevant formulae are 

a{s) = J dxi dx2 $(xi)$(x2) o"o(s) , 

r(i + i/3) 

1 



*(^'?') = ^^TT^h^^^i^ - - + ^) (10) 



+ -f3' 



ln( X ) 

-4(1 + x) ln(l -x) + 3(1 + x)lnx- - 5 - x 

1 — X 



with j3 defined in eq. (^). In order to allow particular QED studies, the 
generation of photons in NEXTCALIBUR has been kept independent for each of 
the two incoming legs. This implies the possibility of two different scales 
and for the electron and positron leg respectively. Mixed situations are 
also possible, in which one leg produces a collinear photon, while the other 
generates a photon with pf. This can be controlled at the level of the input 
file, as explained in the next sections. 

A second issue, relevant when studying QED radiation, is the choice of the 
to be inserted in the Structure Functions of eq. and/or eq. (^^. The 
choice g^ ~ s is proven to reproduce accurately the inclusive four-fermion 
cross sections, at least for s-channel dominated processes. For t-channel 
dominated processes a different choice is more adequate, as studies of certain 
processes have shown. When an exact first order QED radiative correction 
calculation exists for a t-channel dominated process, then the result can 
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be compared to a Structure Function calculation with a scale related to 
the virtuality of the exchanged t-channel photon. With such a value 
the two kinds of calculations agree for small angle Bhabha scattering [ITU 



and multi-peripheral two-photon processes [11|, where the exact calculation 



already exists for some time |jT2|]. When no exact first order QED correction 
calculation is available, the first order soft correction may also serve as a 
guideline to determine [0, 13|. In NEXTCALIBUR, the choice of the scale 
is performed automatically, event by event, according to the selected final 
state, as shown in Table 1. 



Final State 




ql 


No e± 


s 


s 


1 e" 


\t-\ 


s 


1 e+ 


s 


\t+\ 


1 e~ and 1 


\t-\ 


\t+\ 


2 e~ and 2 e"*" 


min( t_ ) 


min(|t+|) 



Table 1: The choice of the QED scale in NEXTCALIBUR. are the scales 
of the incoming while t± represent the t-channel invariants obtained by 
combining initial and final state momenta. When two combinations are 
possible, as in the last entry of the table, that one with the minimum value 
of \t\ is chosen, event by event. 

When using Structure Functions with a t-channel scale, one faces ineffi- 
ciencies in the event generation. To illustrate the problem, we start from a 
simple model, namely a single collinear photon generated with a distribution 
given by the leading behaviour of eq. (0), namely 




The relevant integral is 

/ = dx j d{PS)i^{x)\M{t)\^ , (12) 

where d{PS)4 is the four-body phase-space integration element and |M(t)p 
the t dependent four-fermion Born matrix element squared. The usual way to 
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cure the peaking behaviour of ^(x) would be performing a change of variable 



so that 



[1 — x) = , < p < 1 uniformly (13) 



l^'dpJd{PS),(^-pj (l-a;)(^-^«)|M(t)p, (14) 



and choosing Zg = z, namely 

= - fln(^) - l) (15) 

with tg = \t\. 

The problem here is that, with our QED model, x has to be generated, 
through eqs. (0) and (0), before knowing t. The latter is in fact related 
to the subsequent generation of the four-fermion event. An initial tg must 
therefore be employed, that is not directly related to the true generated t. 
Therefore, event by event, the numerical value of the term 

(l-x)("'^«) = (l-x)^'°^ (16) 

in eq. ([T^) may vary considerably, leading to Monte Carlo inefficiencies. 
Our strategy is as follows. Since eq. ( p^lf ) does not depend on tg, we are 
free to integrate over an extra uniform variable < po < 1, that we also use 
to generate tg distributed as 1/tg between a minimum (t~) and a maximum 
value (t^ = s) 



t 



po 



^. = V(^1 • (17) 
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Our integral then becomes 



/ = f^dp^f^dpjd{PS),[-^'^{l-xt~'^^\M{t)\\ (18) 

with Zg and tg given in eqs. (|^) and (|l^). We have chosen tg distributed 
as 1/tg because this is the expected behaviour of the true t coming from 
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|M(t)p, when small t-cliannel scales dominate. However, this alone does not 
yet prevent a blowing of the term in eq. (jl^) for a few events. We therefore 
split the integration range of po into n equidistant bins, introduce a set of 
weights ttj, normalized such that Y17=i^i — 1; rewrite 



n 



« 1 

i=l 



pdpof^dp J d{PS),a, (^-^) {l-xt''^~^\M{t)\ 



2 



Po. = Po. = -- (19) 

Since Po« < Po < Poi corresponds to the i^^ generation interval 

t- <tg< t+ , t% = t- l^pj , (20) 



eq. (|T9D allows to weight differently events produced in different intervals. 
Notice that the set of weights {ai} can depend on t. A suitable choice is 



-In 7^ 



ai = ai{\tj\ ) = —. (21) 



In 

The meaning of \tj\~ in the above formula is as follows. Once tg is produced 
as described, the four-fermion event can be generated with a given value of 
t (the true t). The variable |t| can then be binned as done for tg in eq. (pO|) 
and, if t'j < \t\ < i^-, \tj\~ = tgj. Then, even by event, the set {a^} changes, 
compensating the term in eq. ([16| ) and improving the generation efficiency. 
Summarizing, the Monte Carlo implementation of the technique is 

1) divide the interval [0,1] in n equidistant bins; 

2) generate po uniformly in [0,1] 
=^ i and t~ are known; 

3) generate ISR and the four-fermion event 

=^ the event weight w, the variable t, j and \tj\~ are known; 
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4) w ^ w ■ n ■ ai{\tj\ ) . 

The algorithm works unchanged when both initial state particles radiate, 
and when photons are generated with a non vanishing pf In this latter case, 
eq. (|l2) is replaced by 

$(a;,c)~- ^-7-^ ^ ^. (22) 

^ ' (l-a;)(i--) i_c + 2^ ^ ' 

The generation of x is as described previously. In addition, c is generated 
with distribution 

1 



l-c+2^ 



and the weights modified accordingly, in order to compensate the term 

l-c+2^ 

% (23) 

appearing in the event weight. 

A last comment is in order on our QED modelling. By using = \t\ 
in the Structure Functions, a few events are produced with very low values 
of ln(|t I /rrig), when \t\ ~ m^. This means that logarithms are no longer 
leading, breaking down the Structure Function approach. Constant terms 
should in principle be included. On the other hand, low values of ln(|t|/mg) 
imply $(x) ~ 6{1 — x), resulting in a suppression of the QED radiation, 
as confirmed by direct calculations [0, |12]. Notice that even implementing 
the exact form of the Structure Functions, also valid for \t\ < ml only 



includes the factorizable part of the missing constants and not all of them. 
In the present version of NEXTCALIBUR we follow a different strategy and 
introduce a minimum |t|mm, such that, if \t\ < \t\min, the ISR logarithm 
is always evaluated at = \t\min f\ With |t|m.m = 100 m^, the minimum 
value of the logarithm still gives physical results, but is small enough to 
reproduce the discussed radiation suppression. We checked the insensitivity 
of our results to the choice of \t\min in the range 50 ml < \t\min < 200 m^. 



We consequently put t = \t\min in eq. (17) 
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5 The running of aEM 



When studying high energy processes, part of the higher order corrections can 
be reabsorbed in the Born approximation by using the so-called Gp scheme. 
In such a scheme Gp, Mz and are input parameters, while the weak 
mixing angle and a em are derived quantities: 

sir = 1-M^/Ml, 

= V2^^^^. (24) 

TT 

In the presence of low t-channel scales such an approach fails, since the choice 
a(t ~ 0) ~ 1/137 is certainly more appropriate for certain sets of diagrams 
and/or kinematical regions. The question is therefore how to include con- 
sistently the running of apM in four-fermion processes with one or more 
electrons in the final state. 

An exact and field-theoretically consistent solution to this problem is rep- 



resented by the Fermion-Loop approach of refs. [|l^]-|]T3) where the whole set 
of fermion one-loop corrections is taken into account by computing running 
couplings g{s) and e(s) and re-summed bosonic propagators. 

In the presence of the WW'j vertex, also loop mediated vertices are re- 
quired to preserve gauge invariance. On the contrary, when the WW'j cou- 
pling is absent, the neutral gauge boson vertices, induced by the fermion loop 



contributions, are separately gauge invariant [0. 

Since the exact expression of the loop-induced vertex functions is rather 
cumbersome, we follow a simplified approach (Modified Fermion-Loop ap- 
proach) by neglecting the separately gauge invariant neutral boson vertices 
and including only the part of the WW'j loop function necessary to pre- 
serve the U{1) gauge invariance. Besides running couplings, we use bosonic 
propagators 



P!i''is) = {s-M'^{s)y' 

pr{s) = {s-MUs)y' 




(25) 
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with running boson masses defined as 



Tw,z{s) are contributions due to the top quark, fi^i^z the complex poles of the 
propagators and 

= St ' = 1 - ^^(^) ■ 

The explicit form of the running functions e^(s), g^{s) and T^y^(s) can be 
found in appendix together with more details on our procedure. Here we 
just mention that, since the leading contributions are in the real part of the 
running couplings, we replace complex poles by complex masses and only 
consider the real part of the corrections. This in practice means replacing 

clil^z) cj{Ml), TwAf^n.,z) ^ TwAM^,z), (27) 

in the above formulae. Note that, as for the treatment of the gauge boson 



widths, this is equivalent to the fixed width scheme [|17 . 

When also the WW'-f coupling contributes, we introduce, in addition, the 
following effective three gauge boson vertex 



p ' 
7m -wvn*^ = i e{s)V^^p 




with s = , s"*" = , s = p^ and 

yi^vp = 9,iu{p-p+)p + 9up{p+-P-),,{l + Sv)+9p^i{P--p)u 
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+ 



{p+ - P- 



s — 



9[s_ 



1 P+uP+p 



g{s+)g{s-){s- - s _ 
+ W^)-g{s-)\\s-g{s^)-,s^g{s-)\ 




(28) 



It is the easy to see that, with the above choice for V^jyp, the gauge 
invariance - namely current conservation - is preserved, even in presence of 
complex masses and running couplings, also with massive final state fermions. 

To preserve U (1), one can either compute g{s) at a fixed scale (for example 
always at s = M^), while keeping only the running of e(s), or let all the 
couplings run at the proper scale. With the first choice the modification 
of the three gauge boson vertex is kept minimal (but the leading running 
effects are still included). With the second choice everything runs, but a 
heavier modification of the Feynman rules is required. Since our approach 
is an effective one, its quality can be judged only by comparing with the 

m 



exact calculation of ref. 



We found that the second choice gives a better 
agreement for leptonic single- final states, while the first one is closer to the 
exact result in the hadronic case, which is phenomenologically more relevant. 
Therefore, we adopted this first option as our default implementation in 
NEXTCALIBUR. 

When using NEXTCALIBUR in the described running coupling mode, U{1) 
gauge invariance is preserved but SU (2) is, in general, violated. The effects of 
such a violation are small at LEP2 energies. In fact our method turned out to 
be numerically equivalent, at -y/i ~ 200 GeV, to the IFL„ approach of ref. [|TS 



's ~ 200 GeV, to the IFL^ approach of ref. 
The authors of ref. [0 divide all Feynman diagrams in gauge invariant sets 
and use different electromagnetic couplings for each set, therefore preserving 
SU(2) X U{1) gauge invariance. Our method is more suitable in all cases 
when no easy separation in different classes of Feynman diagrams is possible. 

An alternative SU{2) x U{1) preserving scheme is represented by the 
formalism of ref. |T^, where extra terms are introduced in the Lagrangian to 
compensate the self-energies contributions and restore gauge invariance. 

A version of NEXTCALIBUR implementing the equations of ref 
rently under development. 
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6 Structure of the program 



We shall briefly describe here the general structure of the program. In the 
MAIN of NEXTCALIBUR the input file is read and various initializations are 
performed, for the matrix element evaluation and the phase space generation. 

The phase space is initialized by calling SUBROUTINE SETPRO, where the 
kinematical channels for the phase space generation are build up, according 
to the chosen four-fermion final state. In SUBROUTINE SETPRO also the input 
parameter set has to be specified. The values to be provided are Gp (GFERMI), 
Mz (ZM), Mw (WM), Vz (ZW), sin^ Ow (SINW2), aEu (ALPHA), and Vw (WW). The 
recommended choice, implemented by default, is the Gp scheme described 
in section ^. A value for the Higgs mass (HM) has to be specified as well and 
the corresponding tree level value of Th (HW) is computed with the formula 



cJ 



rrir = 1.777 GeV, = 0.75 GeV, mi, = 2.9 GeV. (29) 

Finally, SUBROUTINES PHYSICS and HELAC_INIT are called, to initiahze the 
evaluation of the Matrix Element. 

The generation of the QED Initial State Radiation is performed in the 
MAIN of the program. The used algorithm has been extensively described in 



section ^. The number of bins used in eq. ([T9|) to integrate over t-dependent 
Structure Functions is set by the variables nisrlO, nisrll, nisr20 and 
iiisr21. All remaining subroutines are devoted either to the matrix element 
evaluation or to the phase space integration. 

Additional cuts, besides those ones specified in the input file, must be 
implemented directly in SUBROUTINE CUTS, where the two commons 

COMMON/AREAIO/PMI (0 : 4 , 1 : 6) , PM4 (12 : 65) , OMCTl (1:6,3:6) 
COMMON/PHOTONS/PMIG (0 : 4 , 1 : 3) 

contain the fermion four-momenta computed in the Lab frame (PMl), the 
invariant mass squared among all possible particle pairs (PM4), the quantities 
1 - cos% (OMCTl) and the momenta of the two ISR photons (PMIG). If the 
event is rejected LNOT= 1, and the weight is put to zero. 
The conventions for the momenta are as follows: 
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PMl (0 , j ) is the energy of the j*^ particle 

PMl ( 1 , j ) is the X component 

PMl (2 , j ) is the y component 

PMl (3 , j ) is the z component 

PMl (4 , j ) is the four-momentum squared. 



Particles number 1 and 2 are the incoming and e~, respectively. The 
beam is along the x axis and is along the positive direction. 

PM1G(0:4,1) is the four momentum of the photon generated by the in- 
coming e"*", PM1G(0:4,2) that one coming from the e~ and PM1G(0:4,3) is 
the most energetic between the two. PM1G(0:4, j) is different from zero only 
when ISR is generated with a non-vanishing pt distribution. 

The event weight is W, the final value of which is computed immediately 
before the FORTRAN line 

IF (I.GT.O) CALL INB00K(1 , W, UNIT) . 

The status of the run can be checked with the help of the file monitor, 
in which the output file name and the used Monte Carlo points are printed 
out every 5000 iterations (with negative values during the warming phase). 
In the file backup the intermediate result is stored every 50000 Monte Carlo 
points. 



7 Compiling instructions 

NEXTCALIBUR is written in FORTRAN 90. The main program is nextcalibur . f 

and there are five included files, namely 

declare .h 

declare_dp . h (to run in double precision) 
declare_qp.h (to run in quadruple precision) 
compl_mass .h 
list.h. 

NEXTCALIBUR can be run both in double and quadruple precision. Double 
precision is sufficient to run all processes without electrons (or positrons) 
among the final state particles, or processes with final state electrons (or 
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positrons) with at most one final state electron (or positron) allowed in the 
zero-angle forward region. To run in double precision, typical compiling 
instructions are: 

> cp declare_dp .h > declare. h 

> f90 -0 nextcalibur .f -o nextcalibur .out 

> nextcalibur. out < nextcalibur . in 

where nextcalibur . in is the input file (see next section). 

To run processes with more than one final state electron (or positron) 
allowed in the zero-angle forward region, quadruple precision should be used, 
instead. A typical example is the process e'^e~ — > e~^e~fi~^fi~ without any 
cut. To run in quadruple precision, typical compiling instructions are: 

> cp declare_qp.h > declare. h 

> f90 -0 -double_size 128 nextcalibur. f -o nextcalibur. out 

> nextcalibur. out < nextcalibur . in 

The flag -double_size 128 (that treats double precision real numbers as 
quadruple precision real numbers) is essential. 

8 Input 

In this section we list the meaning of the input parameters to be specified by 
the user: 

PAR (3) (CHARACTER* 2) 

Produced fermion with label 3 (to be chosen among 'EL\ 'NE' , 'MU' , 'NM' , 
'TA' , 'NT' , 'UQ' , 'DQ' , 'CQ' , 'SQ' , 'TQ' , 'BQ'). 

PAR (4) (CHARACTER* 2) 
Produced fermion with label 4. 

PAR (5) (CHARACTER*2) 

Produced fermion with label 5. 

PAR(6) (CHARACTER*2) 
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Produced fermion with label 6. 
NIPT (INTEGER) 

The number of points for the Monte Carlo integration. 
NWARM (INTEGER) 

Number of points for the warming up of the phase space generation but 
not used for the actual computation. For processes with zero angle elec- 
trons in the final state it is recommended to chose a non zero value (e.g. 
NWARM= 10000). 

NOPT (INTEGER) 

Number of points for the phase space optimization but also used for the com- 
putation. 

ISPEPMAX (INTEGER) 

Number of iterations for optimizing the a-priori weights. 

OUTPUTNAME (CHARACTER* 15) 

The name of the output file. 

KREL (INTEGER) 

Selects the kinematical channels to be used for the phase space generation. 
KREL=0 is the recommended value for normal runs. 

LQED (INTEGER) 

It includes (1) or excludes (0) ISR. 
LQ_PT1 (INTEGER) 

Collinear QED radiation from the incoming (0) or generation of ISR pho- 
tons with non zero pt (!)■ 

LQ_PT2 (INTEGER) 

Collinear QED radiation from the incoming e~ (0) or generation of ISR pho- 
tons with non zero pt (1). 

LQ_SC1 (INTEGER) 
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t-channel ISR scale for the incoming e+ off (0) or on (1). When choosing (0), 
ISR is computed with the scale = s. When choosing (1), if the process 
has at least 1 final state e+, the used scale is Q'^ — \t\. When there are 2 e+ 
in the final state, the used scale is Q"^ — min{\t\). 

Lq_SC2 (INTEGER) 

t-channel ISR scale for the incoming as above. 
I WIDTH (INTEGER) 

Option for the treatment of boson masses and running of a em- If IWIDTH=1 
all boson masses arc taken to be complex, also in the couplings. When 
IWIDTH=2 aEM is running and the Modified Fermion-Loop is switched on. 

IHIGGS (INTEGER) 

Higgs diagrams included (1) or not (0). 

REL (REAL*8) 
Electron mass. 

RNE (REAL*8) 

Mass of Ug. 

RMU (REAL*8) 

Muon mass. 

RNM (REAL*8) 
Mass of u^. 

RTA (REAL*8) 

Tau mass. 

RNT (REAL*8) 
Mass of Ur- 

RUQ (REAL*8) 
Mass of the u quark. 
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RDQ (REAL*8) 
Mass of the d quark. 

RCQ (REAL*8) 
Mass of the c quark. 

RSQ (REAL*8) 
Mass of the s quark. 

RTQ (REAL*8) 
Mass of the t quark. 

RBQ (REAL*8) 
Mass of the b quark. 

ROOTS (REAL*8) 

The total energy of the colhding e+ and e' . All energies are in GeV. 
SHCUT (REAL*8) 

Minimum value of the invariant mass squared of the event after ISR. 
ECUT(3) (REAL*8) 

Minimum energy of particle number 3. 
ECUT(4) (REAL*8) 

Minimum energy of particle number 4. 
ECUT(5) (REAL*8) 

Minimum energy of particle number 5. 
ECUT(6) (REAL*8) 

Minimum energy of particle number 6. 
SCUT (3, 4) (REAL*8) 

Minimum value of (p(3) + p(4) )^. All invariant masses in GeV^. 
SCUT (3, 5) (REAL*8) 
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Minimum value of (p(3) 



+ p(5))2. 



SCUT (3, 6) (REAL*8) 

Minimum value of ( p(3) 

SCUT (4, 5) (REAL*8) 
Minimum value of (p(4) 

SCUT (4, 6) (REAL*8) 
Minimum value of (p(4) 

SCUT (5, 6) (REAL*8) 
Minimum value of (p(5) 

CMAX(1,3) (REAL*8) 
Maximum value of cos 9 between particle 1 and 3. 

CMAX(1,4) (REAL*8) 

Maximum value of cos 9 between particle 1 and 4. 
CMAX(1,5) (REAL*8) 

Maximum value of cos 9 between particle 1 and 5. 
CMAX(1,6) (REAL*8) 

Maximum value of cos 9 between particle 1 and 6. 
CMAX(2,3) (REAL*8) 

Maximum value of cos 9 between particle 2 and 3. 
CMAX(2,4) (REAL*8) 

Maximum value of cos 9 between particle 2 and 4. 
CMAX(2,5) (REAL*8) 

Maximum value of cos 9 between particle 2 and 5. 
CMAX(2,6) (REAL*8) 



+ P(6) f. 
+ P(6) f. 

+ p(6) r. 
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Maximum value of cos 9 between particle 2 and 6. 
CMAX(3,4) (REAL*8) 

Maximum value of cos 9 between particle 3 and 4. 
CMAX(3,5) (REAL*8) 

Maximum value of cos 9 between particle 3 and 5. 
CMAX(3,6) (REAL*8) 

Maximum value of cos 9 between particle 3 and 6. 
CMAX(4,5) (REAL*8) 

Maximum value of cos 9 between particle 4 and 5. 
CMAX(4,6) (REAL*8) 

Maximum value of cos 9 between particle 4 and 6. 
CMAX(5,6) (REAL*8) 

Maximum value of cos 9 between particle 5 and 6. 

A last remark is in order. When running with the options LQ_PT1= 1 and 
LQ_PT2= 1 for processes with more than one (or e^) in the final state, 
even by cutting out the beam cone the generation efficiency of NEXTCALIBUR 
may be bad, due to genuine 77 like events kicked out from the beam cone by 
large pt photons. If one is interested just in the total cross section, running 
with LQ_PT1 = LQ_PT2= solves the problem, otherwise a reasonable cut on 
the pt of the generated photons must be applied to restore the efficiency. 

9 Test Run Output 

We conclude our description with an example of a calculation that can be 
performed with NEXTCALIBUR. Notice that no cut is present on the final state 
particles. One should be able to reproduce this output within the estimated 
Monte Carlo error. Using an input file as follows 
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el ! par(3) - produced fermion 

ne ! par (4) - produced antifermion 

uq ! par (5) - produced fermion 

dq ! par (6) - produced antifermion 

250000 ! nipt - number of Monte Carlo points 

10000 ! nwarm - points for the warming 

40000 ! nopt - points for a. p. weights optim 

2 ! istepmax - iterations for a. p. weights optim 

output ! outputname - program name 

! krel (0:5) - Channels for MC Mapping 

1 ! Iqed (0,1) - QED off or on 

! lq_ptl (0,1) - finite pt for incoming e(+) 

! Iq_pt2 (0,1) - finite pt for incoming e(-) 

1 ! lq_scl (0,1) - t-chan scale for incoming e(+) 

1 ! Iq_sc2 (0,1) - t-chan scale for incoming e(-) 

2 ! iwidth (1:2) - option for boson masses and running 
! ihiggs (0:1) - Higgs included or not 
0.51099906d-3 ! rel - fermion masses: 

O.dO ! me 

105.658389d-3 ! rmu 

. dO ! rnm 

1.77705d0 ! rta 

O.dO ! mt - be aware that when 

O.OOSdO ! ruq - fermion masses are zero 

O.OlOdO ! rdq - the computation gets faster! 

1.55d0 ! rcq 

150. d-3 ! rsq 

175. dO ! rtq 

4.7d0 ! rbq 

200. dO ! roots - center of mass energy (GeV) 

O.dO ! shcut - cut on inv. mass squared after QED 

O.dO ! ecut_3 - Energy cut for particle 3 

O.dO ! ecut_4 - Energy cut for particle 4 

O.dO ! ecut_5 - Energy cut for particle 5 

O.dO ! ecut_6 - Energy cut for particle 6 

O.dO ! scut_34 - inv. mass cut for the system (34) 

O.dO ! scut_35 - inv. mass cut for the system (35) 
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O.dO ! scut_36 - inv. mass cut for the system (36) 

O.dO ! scut_45 - inv. mass cut for the system (45) 

O.dO ! scut_46 - inv. mass cut for the system (46) 

0. dO ! scut_56 - inv. mass cut for the system (56) 

1. dO ! cmax_13 - cut on cos_(l,3) 
l.dO ! cmax_14 - cut on cos_(l,4) 
l.dO ! cmax_15 - cut on cos_(l,5) 
l.dO ! cmajc_16 - cut on cos_(l,6) 
l.dO ! cmax_23 - cut on cos_(2,3) 
l.dO ! cmax_24 - cut on cos_(2,4) 
l.dO ! cmax_25 - cut on cos_(2,5) 
l.dO ! cmax_26 - cut on cos_(2,6) 
l.dO ! cmax_34 - cut on cos_(3,4) 
l.dO ! cma:x_35 - cut on cos_(3,5) 
l.dO ! cmax_36 - cut on cos_(3,6) 
l.dO ! cmax_45 - cut on cos_(4,5) 
l.dO ! cmax_46 - cut on cos_(4,6) 
l.dO ! cma:x_56 - cut on cos_(5,6) 

and the input values 

Gf = 1.6637 10-^ GeV-^ Mz = 91.1867 GeV 

Mw = 80.35 GeV Vz = 2.4939 GeV 

^ w - ttVs sm dw - 1 - 

we get the following output file 

output 

process : antiel(l) el(2) > el(3) antine(4) uq(5) aiitidq(6) 

This run is with: 

nipt = 250000 

nwarm = 10000 
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nopt 






4UUUU 


istepmax= 




2 











Iqed 






1 


n 1-1-1- i 

lq_px i 









iq_p Lz 









lq_scl 






■A 

1 


±q_scz 






■A 

1 


iwidth 






2 


ihiggs 









rel 




. 


, 510999D-03 


rne 




. 


. OOOOOOD+00 


rmu 




, 


, 105d5oD+00 


rnm 




0, 


, OOOOOOD+00 






, 


, 177705D+01 


rnt 




. 


A/AAAAAT\ I AA 

. OOOOOOD+00 


ruq 




. 


[~AAAAAT\ /~\ O 

. 500000D-02 


rdq 




, 


^ AAAAAr\ /\ H 

, lOOOOOD-Ol 


rcq 




, 


. 155000D+01 


rsq 




. 


. 150000D+00 


rtq 




. 


, 175000D+03 


rbq 




. 


. 470000D+01 


roots 




, 


, 200000D+03 


shcut 




, 


r\r\r\r\f\r\T\ i /^rv 

, OOOOOOD+00 


ecut_3 




. 


, OOOOOOD+00 


ecut_4 




. 


. OOOOOOD+00 


ecut_5 




, 


. OOOOOOD+00 


ecut_6 




0, 


, OOOOOOD+00 


scut_34 




0, 


, OOOOOOD+00 


scut_35 




0. 


.OOOOOOD+00 


scut_36 




0. 


.OOOOOOD+00 


scut_45 




0. 


.OOOOOOD+00 


scut_46 




0, 


, OOOOOOD+00 


scut_56 




0, 


, OOOOOOD+00 


cmax_13 




0. 


. lOOOOOD+01 


cmax_14 




0. 


. lOOOOOD+01 


cmax_15 




0. 


. lOOOOOD+01 



cmax_16 







lOOOOOD+Ol 


cmax_23 


— 





lOOOOOD+01 


cmcLX_24 







lOOOOOD+01 


cincLX_25 







lOOOOOD+01 


cmax_26 







lOOOOOD+01 


cmax_34 







lOOOOOD+01 


cmax_35 







lOOOOOD+01 


cina:x_36 







lOOOOOD+01 


cma:x_45 







lOOOOOD+01 


cmax_46 







lOOOOOD+01 


cmax_56 







lOOOOOD+01 


Input parameters : 


gf ermi 







116637D-04 


zm 







911867D+02 


wm 







803500D+02 


zw 







249390D+01 


WW 







204277D+01 


sin2w 







223558D+00 


1/alpha 







131958D+03 


hm 







120000D+03 


hw 







236865D-02 



permutation 

1 2 3 4 5 6 

1 2 4 3 5 6 

1 2 5 6 3 4 

1 2 6 5 3 4 

1 2 3 4 5 6 

2 1 3 4 5 6 
1 2 5 6 3 4 



Kinematical Diagrams: 

channel 

1: annihi2(wm) 

2: annihi2(wm) 

3: annihi2(wm) 

4: annihi2(wm) 

5: bremf4(wm) 

6 : bremb2 (wm) 

7: conver3(wm,wm) 
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8 


mult i 2 


2 


1 


3 


5 


6 


4 


9 


multiS 


1 


2 


3 


5 


6 


4 


10 


nonabl (wm) 


1 


2 


3 


4 


5 


6 


11 


nonab2 (wm) 


2 


1 


3 


5 


6 


4 


12 


nonab4(wm) 


1 


2 


3 


5 


6 


4 


13 


rambo4 


1 


2 


3 


4 


5 


6 



Differences in the computation 
of the a-priori weights: 



diff( 1 )= 

diff( 2 )= 

diff( 3 )= 



5.38443233306179 
6 . 18288649090395 
2.25190931218116 



a-priori weights: 



1 





301826D-01 


2 





202852D-01 


3 





210927D-01 


4 





327752D-01 


5 





251601D-01 


6 





873505D-01 


7 





299054D+00 


8 





185762D+00 


9 





199223D-02 


10 





162069D+00 


11 





104346D+00 


12 





251601D-01 


13 





477138D-02 



Cross Section (pb) : 
1 0.79733D+00 +/- 0.46D-02 
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A The kinematical channel MULTIl 

We refer to fig. 1 and split the integration over the massive four-body phase 
space into a 3-body decay times a 2-body decay 

J d{PS)i = J rfs45 13(^3,^6,^45)^2(^4,^5) , (30) 

with 

h{P3,P6,Pi5) = J D{p3,ml) D{pQ,ml) D{pi5, 5^{pi2 - Ps -pe -P45) , 
hiP4,P5) = J D{p4,ml) D{p:s,ml) 5'^{pi5 - Pi - P5) , 

D{pi, m- ) = d% 5{p1 - ml) e{Ei) , pij = Pi + pj , S45 = P45 ■ (31) 
The three-body phase space can be rewritten as 

h{P3,P6,P45) = I dx j dy dn-;idVtQ8[F{x,y)], (32) 

where x and y are the reduced energies of the particles 3 and 6 

i?3 = ^x, E, = ^y, 

dfli = d(f)i d cos 6i is the integration element over the solid angle of particle i 
and C36 the cosine of the angle between particles 3 and 6. Finally, 

xy 



F{x, y) = 2 G(x, y)/^H{x, y) - C36 , G{x, y) = M-x-y + -^ 

H{x, y) = (x^ - 4/i3)(?/^ - Afxe) , M = l + /i3 + /i6-/i, 
fn^ s 

/ii = — , n=—, x^ = {1 + 123) - {^+ ./Jl^y . (33) 
s s 



Eq. (|3^) is equivalent to 

h{P3,P6,P45) = 



4 r~ dx I dy dQs dQe , [S{y - y+) + Kv " V-)] , (34) 

ol J2^/JI^ J \r \X^ yj\ 
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where 



oy a 



2/ 



.36V- -4/X3), /3 = 4(l-|)(M- 
7 = 4(M - x)2 + 44 neix^ - 4//3) . 
As for the two-body decay, the following expression holds 

A(a,6,c) = Va^ + + - 2ab - 2ac - 2bc . 
The first four steps of the generation algorithm are as follows: 

1) generate S45 with distribution 

1 

2) generate x with distribution 

1 



1 — X 

3) generate ^3 uniformly and C3 = cos $3 with distribution 



2 



1 E1E3 — TTlg 

03 — 



Os-Cs |Pl||P3 

4) generate 06 uniformly and ce = cos 9e with distribution 



1 



ae-ce 1^21^61 ^ 

where y is the solution of the equation 

2G(xJ)//ff(^+l = 0. 
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The reason for the above choice of ag is that, at this stage of generation, 
Eq is unknown, but its value when cse = —1 (namely when the multi- 
peripheral singularity is more pronounced) can be computed using the con- 
dition F{x,y) = 0, as read from eq. (p^) with c^q = —1. 

From the generated values of C3, Cq, 03 and 06 one then computes 

C36 = cos 6*3 cos 9q + sin 9^ sin 9q cos(03 — 06) , 

to be compared with two possible values (030) obtained from the conditions 
y = in eq. (|3^). There are four possibilities. If C36 7^ c^q the event has 
vanishing weight. If C36 = cj^g 7^ ^36 5 then y = and only the first term in 
eq. (0) contributes. If C36 = C36 7^ c^q, then the second term in eq. (p^) 
has to be chosen. Finally, if C36 = C36 = cj^6) the two terms are picked up 
randomly and the total event weight multiplied by 2. 

The last part of the generation accounts for the t-channel exchanged 
fermion in MULTIl. By defining q = Pi — Ps, the variable C4 = cos ^4 peaks, 
in the rest frame of S45, as follows 

1 ^ 2go^4 - ^ g45 + ml-ml ^^^^ 



a — Ci' 2\q\\pi\ ' 2y^S45 

On the other hand, also the diagram obtained when replacing p^ <-> p^ con- 
tributes with an analogous peaking behaviour 

1 ^ _ 2go^5 -q^ S45 + ml-ml ^^^^ 



Since, in the rest frame of S45, C4 = — C5, both distributions can be simulta- 
neously mapped, via the last three items of our generation algorithm: 

5) go the rest frame of S45; 

6) generate 04 uniformly and C4 with distribution 

1 



(a - C4){b + C4) 
7) boost p4 and p^ back to the lab frame. 
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A massive four-fermion event is then generated, taking into account all peaks 
due to t-channel exchanged massless or nearly massless particles. 

In NEXTCALIBUR the described procedure is further refined by making it 
fully symmetric under the replacements {^1,^3} ^ {P2,P6}- 

Finally, the same ingredients are used also when constructing the chan- 
nel MULTI2, the only variation being a different choice of the distributions 
employed in the generation of some of the variables. 



B The Modified Fermion-Loop approach 

Following ref. [|T^ we split the bare W and Z boson self-energies into universal 
and non universal contributions 

-2 



^z{p' 



tzip') + %Tz{p 



(41) 



where g is the bare weak coupling constant s the sine of the bare mixing angle 
and = 1 — . When including fermion loops only, the explicit expressions 
for the above quantities are 

f 1 



487r2 2 



9^ 
48c%2 



Y.N'f[l{p^) + {l + 2hfYf)F{p\mf)} , (42) 
/ 

E^/{(^? + 4)A/) 



+ {v) + a] + AhfYfa])F{p\mf)} 



with Yf = 2{Qf — I^f), Vf = I-if — 2s^Qf , aj = —/a/ and Nj = number of 
colors. Furthermore 



Hp') 

F{p'^, m) 
G{s, t, m) 

m 



- + ln( V) - ^ 
e 6 



-p'fie 
-6pV(e) G{p^,p^,m 



I f - sx{l - x) 

ax { x[l — X) m ; : — 

I ^ ^ -txil-x) 



vraF 1 - - 
2 



with e = — 4 . 



(43) 
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Finally, the leading contribution to Tw^zip"^) originate from the top quark: 



Twip') = /(e) 



167r2 



m. 



- 2p^ 



Mb) - Mb) + 4 - iln(6) - |G(p^p^m,) 

-30 D 6 



167r2 



-m 



- + Inm^ + Jo(6+) + Mb-) 



+ p^ G{p^,p^,mt)^ 



(44) 



where 



J.(6) = dx x' ln(l + bx), b= , 6± = — , xl-x± + l= 0. (45) 

Two more expressions are necessary, namely the photon (S-^) and the photon- 
Z mixing (Sx) self-energies 



127r2 



j:Q}N]{np') + F{p^mf)} , 



247r2 c 



^E^/Q/^/U(p') + ^(p'-^/)}- 



(46) 
(47) 



In terms of self-energies the running couplings are defined as follows 



, Mp') 



9^{p'^) ' 



(48) 



Since, as explained in section ^ we are only interested in the real part of the 
corrections, we renormalize the W and Z boson propagators at the real mass 



Ml 



fly, - ReSi^(M^) , 
/i,-ReZ(M^), 



where fi^^z are the bare gauge boson masses squared and 

tW) 



Z{p') = S^(p^ 



p2 + S^(p2 



(49) 



(50) 
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Finally, Gp fixes our third fitting equation 



V2 Sifiw-f^wm' 



(51) 



By solving the above equations it is easy to get g^, fi^ and in terms of 
and M|. Putting everything together, the running couphngs read 



g^s) = M^aRe \ l + -M^—G{M^,s,0)- 



and 



eHs) 



+ 



where 



[G{s,s,me) 



^'i- ^-yy^-^-' 87^2 

+ G{s, s, m^) + G{s, s, rrir) + 2 G{s, s, my) + 2 G{s^ s, rric) 
+ 2 G{s, s, mt) - 2 ^(M^, M^, m^) + ^(s, s, m<i) 

+ ^(s, s, rus) + ^(s, s, mb)] - a [Tw{M^) - TV(0)] , (52) 



M^aNRe 



Q /\/f2 ^ 

1 - a[TH.(M^) - Twm + ^^G(M^, M|, 0) 



2 7r2 

2 



(G(M^, M^, m,) - G(M|, M|, m,)) 



TT 



1 



G(s, M|, 0) + -[G(s, s, me) + ^(s, s, m^) 



1, 



+ ^(s, s, rrir)] + g [^^(-5, -s, m„) + G(s, s, mc) + G(s, s, m*)] 

- ^G(M|, M|, mt) + ^^[^(s, s, ma) + ^(s, s, m,)] 
D 24 



p. 



(53) 



N 



1 - 



Ml 



l-aRe(Tz(M|)-Ti^(0)) 



a Re (T;y(M^) - T^(0)) + ^^G(M^, M|, 0) 
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+ 



4 7r2 



Re 



G{M^, M^, mt) - G{Ml, Ml, nit] 



(54) 



A comment is in order on the appearance of the hght quark masses in the 
above expressions. The weak couphng g'^{s) gives sizable contributions at 
high energies, where all light fermion masses can be neglected, so that only 
the top terms survive in eq. On the contrary, being interested in the 

low energy regime of eq. (|53|) , we have to include in our formulae the hadron 
vacuum polarization A^(s), as computed in ref. [^, rather than using quark 
masses. By inserting 



m eq. 

e\s) 



Re G{s, s, m) 

J{s, m) 
, one gets 



M^aNRe 



I 5 „ N 

-In — + — + Re J(s,m) 
\s\ 18 



dxxil — x) In ( 1 



:x{l 



X 



(55) 



a[Tw{M^) - Twm + l^G{M^, M|, 0) 



+ 



+ 



\^ (G(M^, M^, m,) - G(M|, Ml, 
2 



2 7r2 

2 



37r2 



M^aN J{Ml,mt)- J{s,mt) 



25 1 mj I ^ 

— > In — I- H > In 

^« ^ Ml 36 ^ 



108 48 



+ 



1 „ 1 



q=u,c 



IT 



M| 



144 



g=d,s,i 



i=e,iJ,,T 



16«(0) 



(56) 



with l/a{0) = 137.0359895. Notice that, due to our choice of the input 
parameter set, a dependence on the light quark masses is left in the previous 
equation. However, since the pure bosonic corrections are missing, mt can 
be interpreted as a free parameter |T6[ to be fitted in such a way that a(0) = 
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e^(0)/47r = 1/137.0359895. The quark masses appearing in eq. (^) only- 
affect this fitting procedure, and not the running of the electro-magnetic 
coupling. In other words, a change in the m^'s is compensated by a different 
fitted value for m^. For consistency, in NEXTCALIBUR the light quark masses 
are taken to be those that reproduce the hadronic vacuum polarization, but 
the numerical results are rather insensitive to the initial choice of m^'s. 

Eqs. ( ^21) and (^6|), together with equations (p5|)-(p8l) completely define 
our Modified Fermion-Loop approach. 
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